function [P] = Func_Extract_Field(f, distance, TX, RX, shape, cross_section, dx, ds, E)
    % Shap definitions used in Func_VPE
    if shape == 1
        H_top = cross_section{2}/2;
        W_half = cross_section{1}/2;
        H_bottom = H_top;
        
    elseif shape == 2
        H_top = cross_section{2};
        W_half = cross_section{1};
        H_bottom = cross_section{4};
    elseif shape == 3
        H_top = cross_section{2};
        W_half = cross_section{1};
        H_bottom = cross_section{4};
    elseif shape == 4
        H_top = cross_section{2};
        W_half = cross_section{1};
        H_bottom = cross_section{4};
    end

    dy = dx;                        % cellsize along vertical direction
    nh_top = ceil(H_top/dy);        % number of cells along the top part, veritically
    nw_half = ceil(W_half/dx);      % number of cells for half width part, horizontally
    Nx = 2*nw_half;                 % number of cells along the horizontal direction
    Nz = ceil(distance/ds);         % number of cells along the propagation direction
    x_c = Nx/2 + 1;
    y_c = nh_top + 1;

    Rx_horizontal = RX{1};
    Rx_height = RX{2};
    Rx_vertical = Rx_height - H_bottom;
    rp_y = y_c - round(Rx_vertical/dy);
    rp_x = x_c + round(Rx_horizontal/dx);
    E_extract = zeros(1, Nz+1);
    E_extract(1:end) = E(rp_y,rp_x,:);

    c = 299792458.0;
    lambda = c/f;
    mu_0 = 4*pi*1e-7;
    eps_0 = 1/mu_0/c/c;
    eita = 1/sqrt(eps_0/mu_0);
    l_eff = lambda/pi;
    Rs = 73;
    E_dB = 10*log10( (l_eff.*abs(E_extract)).^2/8/Rs );
    Tx_Gain_Diff = (TX{4} - 30) + 10*log(eita/2/pi) - TX{3};
    Rx_Gain_Diff = RX{3} - 2.15;
    VPE_shift = 30;
    P = E_dB + Tx_Gain_Diff + Rx_Gain_Diff + VPE_shift;
    P = P(1:round(RX{4}/ds):end);
end